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ABSTRACT 

A restorative-iterative proceedure is tested in initializing the 
height and wind fields for a global barotropic prediction model. Real 
500 mb data is used as the basis of initial height values in the mid 
latitudes and as the basis of the initial streamfunction values in the 
tropics. The proceedure is tested for the case of frictionless and 
adiabatic flow and for the case where friction and a simulated heating 
function is included. In addition, experiments were performed in which 
an n observed M divergent wind was added in the tropics. Results are very 
encouraging in the sense that the amplitudes of the external inertia- 
gravity motions excited due to initial imbalance between the mass and 
wind fields are reduced considerably with the use of the proceedure. 
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ABSTRACT 



A restorative-iterative proceedure is tested in initializing the 
height and wind fields for a global barotropic prediction model* Real 
500 mb data is used as the basis of initial height values in the mid 
latitudes and as the basis of the initial streamfunction values in the 
tropics. The proceedure is tested for the case of frictionless and 
adiabatic flow and for the case where friction and a simulated heating 
function is included. In addition, experiments were performed in which 
an "observed" divergent wind was added in the tropics. Results are very 
encouraging in the sense that the amplitudes of the external inertia- 
gravity motions excited due to initial imbalance between the mass and 
wind fields are reduced considerably with the use of the proceedure. 
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1. Operational numerical weather prediction is on the threshold of using 
rather sophisticated global prediction models. This is being made feasible 
by the rapid expansion of computing capacity as well as development related 
to general circulation research. But unlike the general circulation model 
the operational model must begin from real data. For the foreseeable 
future the initial data will be somewhat fragmentary in nature so this 
naturally poses some difficult problems. For example, one such problem 
might be; given a certain type of data, i.e., observations of mass alone, 
wind alone, or some mixture, what method should be used to initialize the 
missing type of data so that a global model will not suffer during the early 
part of a forecast run from excitation of excessive inertia-gravity motions? 
This paper will be concerned with this topic. 

As a general guide the geostrophic adjustment theory indicates that for 
mid-latitude motions of synoptic scale the wind does have a strong tendency 
to adjust to a given or observed mass field. In other words strong changes 
takes place in the wind field rather than the mass field as the adjustment 
is occuring. In the tropics rather the opposite may occur. For a theoretical 
treatment of these subjects see (Cahn, 1945, Obukhov, 1949, Bolin, 1953, or 
Washington, 1964). From quasi-geostrophic theory, most valid certainly in 
the higher latitudes, certain balance relationships are derived which allow 
a computation of such a balanced wind if the mass observations are adequate. 
Both non-divergent and divergent parts of the wind can be computed. On the 
other hand a mass field may be computed from a non-divergent wind if it is 
adequately observed. This may be applicable in the lower latitudes. 

Such computed initial winds are presently being used in operational 
mid-latitude models (Shuman and Hovermale, 1968, Kesel and Winninghoff, 1971) 
but are not adequate for the tropics. Also difficulties do arise in the 
computation of such winds. The equations in the most sophisticated technique - 
the total non-linear balance equation and the consistent form of the omega 
equation - are not easy to solve. There is the ellipticity difficulty and 
if another coordinate system other than the pressure system is used in the 
vertical much complication arises. 

Therefore in the past five years another general method, proposals for 
which have been discussed by Nitta and Hovermale, 1969, Miyakoda and Moyer, 
1968, Winninghoff, 1968, Nitta, 1969, Chamey, Halem and Jastrow, 1969, 
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Jastrow and Halem, 1970, and Mesinger, 1971, has been proposed as an 
alternate method of achieving a total balanced condition. This method which 
we will refer to as the "restorative-iterative" method, though all of the 
above references do not use explicit restoration, is one in which the 
principle involved is the use of the equations in an iterative sense either 
at a fixed time level or even in a four dimensional sense in which data is 
assimilated into a running model. This allows for the natural adjustment 
mechanism itself to achieve the desired balance or a time scheme such as 
the Euler-backward scheme may be used for the same purpose. The method has 
the advantage of mathematical simplicity and complete consistency with the 
prediction model though at present it is expensive in computer time. This 
disadvantage, however, will very soon be removed by computers such as the 
ILL I AC IV. 

This paper will be concerned with this restorative-iterative proceedure 
but only at a single time level. However, it will be applied to the entire 
globe. 

2. Model 

The experiments to be reported here used a simple free surface barotropic 
prediction model. The equations for this model are, 
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Equations (1) and (2) are the zonal and meridional momentum equations in 
flux form. Equation (3) is the continuity equation. Symbols are, 
u - zonal component of wind 
v - meridional component of wind 
h - height of the free surface 
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d\ - increment in longitude 
a - radius of the earth 

- increment in latitude 

- latitude 

f - coriolis parameter 
g - acceleration due to gravity. 



The finite difference equations corresponding to equations (1), (2), and 
^3) are written on a staggered grid where the variables are located relative 
to one another as shown in figure 1, with h alone carried at the poles. 

The increment of latitude and longitude used was five degrees. Therefore, 
there are 1260 points (36 x 35) to cover the globe. The differencing used 
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the experiments. The spatial differencing 
is shown in Appendix A. 

The time step used is 15 minutes. 

This was possible in spite of the 40 kilo- 
meter distance between grid points at 
85° N and S, because of a proceedure to 



average the effect of inertia-gravity waves of high frequency in the zonal 
direction used by Arakawa (Langlois and Kwok, 1969). 

A coefficient is defined, 

A. = . 125(D . -1) /D. (4) 
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1 and 2 uh in equation (3) and in the advective terms of equations (1) and 
(2) and h in the pressure gradient term in equation (1) is replaced by, 
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is used. 

3. Experimental Results 

A number of experiments were conducted with the global barotropic model 
to test the effects of a restorative-iterative initialization proceedure. 
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The basic data used was a 500 mb. height field read from the Fleet Numerical 
Weather Central 1 s master file for 00Z, January 8, 1971. This field was 
analized by the FNWC objective scheme. A program was used to .convert it so 
values at 5° latitude-longitude intersections were available over the 
Northern Hemisphere. For the Southern Hemisphere a mirror image of the 
field was used. 

The first result, shown in figure 2, is a plot of the 500 mb. height 
at the poles for a 36 hour forecast in which the wind was initialized 
everywhere by use of. 



\|r = gh/.OOOl 
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v = 
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where i ) \ is the streamfunction and where these computations were done over 
the Northern Hemisphere and reflected into the Southern Hemisphere. The 
wind was set at 0 along the equator. Figure 3 shows the plot of the height 
values at a selected equatorial point. Immediately it is obvious that high 
to moderate frequency oscillations of large amplitude occur in the forecast. 
Realizing that straight use of the leapfrog time scheme is subject to the 
type of instability known as solution separation these runs were rerun 
using the stable two-level Euler-backward every 24 time steps to recycle or 
n tie together" the solutions. Improvement was noted at the poles but not 
enough to avoid the conclusion that the results are completely unacceptable 
for an operational model. The crude initialization here is at fault. 

In figures 4 and 5 we show some similar results for runs in which the 
initialization was improved. The upper curves represent the cases in which 
the linear balance equations. 



V 2 f + V|.Vf/f = gV 2 h/f 
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V 2 h = (fV 2 f + Vf /g 


(11) 



were used to define the initial fields. Equation 10 was used from latitude 
30° to the pole. Equation 11 was used from latitude 20° to the equator. At 
latitude 25° a blend was used. For "observed" \|r values in the tropics equation 
(7) was used. This is permissible as this work is not concerned with veri- 
fication of forecasts per se, but rather the question of initialization. 
Comparison of the upper curves of figures 4 and 5 to the curves in figures 2 
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and 3 shows immediately improvement in the sense that the amplitudes of the 
higher frequency motions are only one-third as great as they were in the 
former runs. 

The lower curves in figures 4 and 5 represent the situation where an 
iterative proceedure was used in addition to the solution of equations (10) 
and (11) to initialize the forecast. What was done here simply was to run 
the equations alternately forward and backward six steps (now properly 
"iteration” rather than "time" steps), allowing the variables to freely 
change. Also during this proceedure the Euler-backward time differencing 
scheme was used. No restoration was used here so no constraint is placed 
on the adjustment mechanism. The process was run the equivalent of only 18 
hours. Note the initial height value for the equatorial point changes about 
seven meters. The results are excellent. 

However, the situations shown in these curves are still rather ideal. 

In order to test the process under more realistic conditions several experi- 
ments were designed to simulate the incorporation of observed data, to 
simulate the effects of large scale heating and friction and to include an 
"observed” wind containing divergence in the tropics. 

The incorporation of observed data gives rise to the restorative feature 
of the iterative process. The results shown in the lower curves of figures 
4 and 5, though excellent, do suffer from the lack of restraint during the 
iterative proceedure since the initial "observed" values of the mass or wind 
may be changed more than desired. Therefore, during the process, after 
each normal iteration of equations (1), (2), and (3), the following is added, 

uh = (1-k )uh + k hu (12) 

u u o 

vh = (1-k )vh + k vh (13) 

V VO 

h = (1-k, )h* + k u h (14) 

h ho 

where the asterisk refers to the values of the variables computed from 
equations (1), (2), and (3) and the subscript o refers to the "observed" 
value of the variable. The k's are a function of latitude, k^ and k^ 

were set at .5 from latitude 20° S to 20° N, 0 from 40° to the pole and a 

linear variation between 0 and .5 between latitude 20° and 40 . k^ was 
set at .3 above and below latitude 40° N and S, 0 between latitude 20° S and 
20° N and a linear variation between 20° and 40°. Figure 6 shows the result 
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of the forecast run. The upper curve is for the polar point(s), the lower 
curve for the same equatorial point as was used in the previous figures. A 
very slight deterioration in the results when compared to the lower curves 
in figures 4 and 5 is apparent. This is expected because restraint is being 
placed on the efficiency of the adjustment. Nevertheless, these results are 
still far superior to the results represented by the upper curves of figures 
4 and 5 . 

Some comment is necessary concerning the particular coefficients used. 
Independent runs by the author and others have indicated that for a baro- 
tropic model in which only the very high frequency external gravity waves 
are present too much interference with the adjustment mechanism is undesirable 
as only the damping effect of the Euler-baclcward time scheme may be retained 
on the high frequency motions. Therefore, coefficients greater than .5 are 
not recommended. But with baroclinic models (see the reference to Mesinger) , 
higher coefficients may be allowable as the dominant internal inertia-gravity 
waves are of lower frequency. 

The next set of experiments concerned the inclusion of heating and 
friction effects into the initialization proceedure. Strictly speaking there 
can be no heating with a barotropic model but the effect may be simulated 
by a systematic addition or subtraction of height as a function of latitude 
or some other factor at each time step. For this we used, 

H = 2 x 10 ^(h) (.5 - |sin (j)|) . (15) 

For friction, 

F = -.00003 u , F = - .00003 v (16) 

u v 

was used. In figure 7 results from a run initialized by use of equations (10) 
and (11) only are shown in which the expressions in (15) and (16) were added 
to the forecast equations. The amplitudes of the high frequency motions are 
very similar to those in the upper curves of figures 4 and 5. It is true 
that here total energy is apparently declining as the polar height rises over 
the 36 hour period and the equatorial height falls but again as the concern 
of these experiments is the reduction of the amplitude of the high frequency 
motions this decline in the total energy is not crucial. In figure 8 results 
are shown where the initialization was done as in the case shown in figure 6. 
No attempt is made to include the effects of the heating or friction terms 
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in the initialization process, that is to initialize the large scale 
divergence which would be in balance with these terms . Except for the 
effect of the heating and friction during the forecast there is nothing 
new here. These curves reproduce essentially the curves from figure 6. 

To include the effects of the heating and friction terms into the 
restorative-iterative proceedure it must be understood that the processes 
of geostrophic adjustment and adjustment to the heating and friction terms 
are somewhat different in nature. For example, the equations can be run 
alternately forward and backward since the natural adjustment mechanism is 
dependent on the dispersive nature of the gravity-inertial waves which 
travel in all directions. Advection is not of prime importance. With 
heating and friction, however, there is no reason to run forward and back- 
ward as this would cancel the effects of these external forces completely. 
Therefore, some splitting of the iterative proceedure is strongly indicated. 
Note also that in the quasi-geostrophic theory the form of the omega equation 
suggests that the effect of heating and friction might be solved for 
completely separate from the part of the process in which geostrophic adjust- 
ment is important. Several variations of such a divided process were tried. 
The best result was obtained when first the restorative-iterative process 
was run as described before, but then before the forecast was begun the 
finite difference versions of. 



d(uh)/dt = - kuh 


(17) 


d(vh)/dt = - kvh 


(18) 


dh/dt + hV-W = H;k = + .00003 


(19) 



were run in a forward sense only with the restoration being done as indicated 
in equations (14), (15), and (16) after each step. k^, k^, and k^ were 
set equal to .5 at all latitudes during this part of the process and the 
restoration was made not to the original values of the variables but to the 
balanced values achieved at the conclusion of the adiabatic part of the 
proceedure. Only three hours of additional iterations were needed before 
no further changes in the height values were noted. Figure 9 shows results 
of the forecast. Differences with figure 8 are quite small though in the 
print-out of the entire field differences near five meters are noticeable at 
many points from the run in which the initialization was completely adiabatic. 
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Finally one question remains. If observations of tropical winds are 
adequate to be meaningfully analysed as has been assumed here, in reality 
they may contain some divergence not in any balance with the mass. That 
is the time change term may well be of the same order of magnitude as are 
other terms in the equations. Quasi-geostrophic theory, in fact, gives us 
little guidance in the low latitudes concerning these matters. 

We created arbitrarily an "observed” divergent wind in the tropics by 
means of, 

uh d = 2h(cos tt(I-1) 19 ) ( co s tt(18-J)/10) (20) 

where I ranges from 1 to 36 and the equatorial value of J is 18. Above 

and below latitudes 25° N and S uh^ was set at 0. In figure 10 the 

results at our points are shown where the model was initialized as for the 

case shown in figure 9, but at the commencement of the forecast uh, was 

a 

added to uh . The curve for the polar point (s) is identical to the polar 
curve in figure 9, but the curve for the equatorial point is considerably 
changed. The amplitude of the large oscillation between hour 12 and 28 is 
nearly doubled. Of course, this result is not unexpected. A question of 
convenience now arises. To what extent, if any, is it necessary to divide the 
observed wind into a divergent and non-divergent part to use the initializa- 
tion proceedure in the tropics? In figure 10 the process was run with the 
restoration made to the non-divergent part of the wind during the adiabatic 
part of the process. In figure 11, on the other hand, results are shown in 
which the initialization proceedure was started with the total wind as 
defined from equations (20) and (7) and in which the restoration used the 
total wind. Several variations were tried as well, but all of the results 
varied negligibly from those shown in figure 11 and figure 10. Therefore, 
we have verified that our restorative-iterative process remains very 
convenient in the sense that total winds may be used without any preprocessing 
into divergent and non-divergent components. 

4. Conclusions 

We conclude from these experiments that a mathematically simple proceedure 
for initializing global primitive equation models is quite feasible given 
adequate observations of the mass in the temperate latitudes and the wind in 
the tropical Latitudes. The proceedure is entirely consistent with the 
forecast equations and can be extended easily to baroclinic models. It can 
be extended to the four dimensional situation. There is no reason why this 
proceedure can not take into consideration a mixed situation in which both 
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wind and mass observations exist in a given area but may not be adequate 
over some entire region. The restoration coefficients can be varied. 

Thus if observed winds were analysed over the U.S. and Canada they could 
be taken into account in that particular area. The whole proceedure will 
soon be feasible operationally as computer technology is so rapidly 
advancing. 
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Appendix A 

Finite Difference Equations 



The difference equations corresponding to equations (1), (2), and 
(3) are. 
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In these equations the following convention has been used* In Fig {A-l) 

the configuration and indexing of the 
mass and wind variables is displayed. 
Therefore, the index X in equations 
(A. 2) and (A. 3) is i itself if j 
is odd and i+l if j is even. 

Also the following definitions apply. 
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Legends for Figures 



Figure 2. Plot of 500 mb height at pole during 36 hour forecast initialized 
everywhere with first order quasi-geostrophic wind and analysed 500 mb 
height field. 

Figure 3. Plot of 500 mb height at selected equatorial point during 36 hour 
forecast initialized everywhere with first order quasi-geostrophic wind and 
analysed 500 mb height field. 

Figure 4. Upper curve, plot of 500 mb height at pole during 36 hour fore- 
cast initialized by solution of linear balance equations. Lower curve, plot 

of same but forecast initialized with addition of iterative proceedure with 
no restoration. 

Figure 5. Same as in figure 4 but for selected equatorial point. 

Figure 6. Upper curve, plot of 500 mb height at pole during 36 hour fore- 

cast in which restoration feature was added to the initialization proceedure. 
Lower curve, plot of 500 mb height at the selected equatorial point for same 
forecast. 

Figure 7. Upper curve, plot of 500 mb height at selected equatorial point 
during 36 hour forecast initialized by solution of linear balance equations 

and in which friction and simulated heating is included. Lower curve, same 
as above for the polar point. 

Figure 8. Upper curve, same as upper curve in figure 7 except initialization 
includes the restorative- iterative proceedure applied adiabatically. Lower 
curve, same as above for the polar point. 

Figure 9. Upper curve, same as upper curve in figure 8 except that the re- 
storative-iterative proceedure includes a diabatic part as described in text. 
Lower curve, same as above for polar point. 

Figure 10. Upper curve, plot of 500 mb height at selected equatorial point 
initialized as in figure 9. "Observed” divergence added to wind in tropics. 
Lower curve, same as above for polar point. 

Figure 11. Upper curve, same as upper curve in figure 14 except the divergent 
part of observed tropical wind was included in restorative-iterative initiali- 
zation proceedure. Lower curve, same as above for polar point. 
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